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ABSTRACT 

The stability of the dynamical trajectories of softened spherical grav- 
itational systems is examined, both in the case of the full TV-body 
problem and that of trajectories moving in the gravitational field of 
non-interacting background particles. In the latter case, for TV > 10000, 
some trajectories, even if unstable, had exceedingly long diffusion times, 
which correlated with the characteristic e-folding timescale of the insta- 
bility. For trajectories of TV « 100000 systems this timescale could be 
arbitrarily large — and thus appear to correspond to regular orbits. 
For centrally concentrated systems, low angular momentum trajecto- 
ries were found to be systematically more unstable. This phenomenon 
is analogous to the well known case of trajectories in generic centrally 
concentrated non-spherical smooth systems, where eccentric trajecto- 
ries are found to be chaotic. The exponentiation times also correlate 
with the conservation of the angular momenta along the trajectories. 
For TV up to a few hundred, the instability timescales of TV-body sys- 
tems and their variation with particle number are similar to those of 
the most chaotic trajectories in inhomogeneous non-interacting systems. 
For larger TV (up to a few thousand) the values of the these timescales 
were found to saturate, increasing significantly more slowly with TV. 
We attribute this to collective effects in the fully self-gravitating prob- 
lem, which are apparent in the time-variations of the time dependent 
Liapunov exponents. The results presented here go some way towards 
resolving the long standing apparent paradoxes concerning the local 
instability of trajectories. This now appears to be a manifestation of 
mechanisms driving evolution in gravitational systems and their inter- 
actions — and may thus be a useful diagnostic of such processes. 

Key words: Celestial mechanics, stellar dynamics - Instabilities - 
Gravitation - Galaxies: kinematic & dynamics 



1 INTRODUCTION 

A persisting and perplexing question regarding the 
dynamics of TV-body systems concerns the interpre- 
tation of the local divergence of their trajectories 
and the origin of this instability. Does this phe- 
nomenon, first pointed out by Miller (1964), tell us 

* Currently at the Center for Astrophysics and Space 
Astronomy, University of Colorado, Boulder 



anything non-trivial about the behaviour of grav- 
itational systems? Does it, for example, have any 
relation with processes, whether collective or local, 
leading to evolution in these systems? Is it useful in 
understanding and characterising such processes? 

A series of studies by Kandrup and collabora- 
tors (e.g., Kandrup et al. 1994) have shown the lo- 
cal instability of trajectories to be an extremely ro- 
bust phenomenon — occurring for all initial condi- 
tions and perturbations investigated, for a variety 
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of initial macroscopic configurations, and with an 
exponential timescale that is roughly independent 
of these properties and of TV, whether the system 
is initially in dynamical equilibrium or not. In po- 
tentials resulting from collections of fixed point par- 
ticles the exponential divergence and its constant 
e-folding time have been checked for TV up to a mil- 
lion particles (Valluri & Merritt 1999). However, as 
might be expected, conservation of the action vari- 
ables, for systems with separable smoothed out po- 
tential, improves with increasing TV. Thus, at least 
in some respects and over sufficiently short time- 
scales, the trajectories of particles in high-TV dis- 
crete systems resemble those in continuous ones — 
since, qualitatively, the motion is completely deter- 
mined by those action variables. This implies that 
time averages of trajectories in steady state discrete 
systems will approximate those of the correspond- 
ing smooth systems for times much larger than the 
exponentiation time-scale (which is a fraction of the 
dynamical time). Nevertheless, the latter do not ex- 
hibit exponential divergence of nearby trajectories 
in the case of separable potentials. Moreover, tra- 
jectories in fixed smooth potentials are characteris- 
tics of the time-independent collisionless Boltzmann 
equation (CBE). If therefore, as is often assumed, 
large TV-body gravitational systems are adequately 
described by the steady state CBE for timescales 
long compared to their dynamical time, the the phe- 
nomenon of the TV-invariant e-folding time must be 
considered paradoxial. 

It can be shown (e.g., Braun & Hepp 1977; 
Spohn 1980) that, for large large-TV gravitational 
systems, the dynamics described by the full equa- 
tions of motion, converge towards that resulting 
from solutions of the CBE, provided that the po- 
tential has bounded first and second derivatives (a 
condition which enters, for example, in the form of 
Eq. 2.52 of Spohn). This amounts to "softening" the 
potential to get rid of the singularity at the origin. 
Indeed, in potentials with a short range cutoff, the 
force due to the mean field is always greater than 
the force between any two particles for large enough 
TV. For example, in a spherical system, this simply 
requires that for a test particle at radius R 



where e is the (Plummer) softening length. Or equiv- 
alently that e > R/yN. This in fact justifies the 
mean field approach embodied in the CBE. 

Because, in the case of point mass systems, the 
binary force between neighbours can be larger than 
the force due to the mean field, formally speaking, 
the dynamics of even a large- TV system is not equiv- 
alent to the one described by CBE — and thus the 
exponential divergence does not in itself signal any 
inconsistency. Perhaps the only puzzling point then 



is why the timescale of the exponential divergence 
does not correlate with particle number while other 
dynamical properties of gravitational systems obvi- 
ously do. (As noted above, this is even true for tra- 
jectories of systems of fixed point particles, where 
the exponential instability cannot be attributed to 
any time-dependence: Valluri & Merritt 1999) . 

Goodman, Heggie & Hut (GHH) conducted an 
extensive investigation of the exponential instabil- 
ity and concluded that it is precisely at the radius 
where the mean field force and the force due to near 
neighbours become comparable that the latter inter- 
actions have the greatest contribution to the process. 
Since the interparticle spacing, projected on two di- 
mensions, decreases as l/VTV, a test particle is likely 
to suffer one such encounter at every crossing of the 
system — independent of TV. For larger TV these en- 
counters last for shorter times and thus the deflec- 
tion angle decreases (roughly as l/\/TV). Neverthe- 
less, the initial divergence between two trajectories 
(in the linear approximation) is independent of TV 
(see GHH). Thus the trajectories of two particles, 
initially separated by a small distance compared 
to the impact parameter, will diverge exponentially 
with a timescale independent of TV. Consider how- 
ever two trajectories separated in such a way that 
one is affected by an encounter with a field star at 
impact parameter ~ i?/\/TV and the other not sig- 
nificantly affected — being at a distance ^> R/yN 
from the perturbing object. The affected test parti- 
cle will be deflected by an ang le 9 ~ l/VTV. Thus, 
for large (compared to R/y/N) separations and large 
TV, the divergence will be small. 

In conventional dynamical systems where chaos 
is present for most initial conditions, it is usually 
either produced by the bulk properties of the po- 
tential or by local large angle scatterings. In both 
cases large deviations of trajectories, resulting from 
the non-integrable component of the potential, are 
expected in times comparable to the exponentiation 
times. In unsoftened gravitational systems on the 
other hand, the contribution to the potential re- 
sponsible for exponential divergence on very short 
timescales independent of TV comes from encounters 
whose range, duration, and deflection angles pro- 
gressively become smaller at larger TV. Their effect 
on the actual trajectories becomes small, and the 
divergence does not lead to any macroscopic evolu- 
tion on the divergence timescale — since it is only 
at small separation that the divergence is indepen- 
dent of TV. Any evolutionary effects due to the expo- 
nential instability, therefore, will have to come from 
larger range interactions, whether individual or col- 
lective. 

In order to answer the questions posed in the 
opening paragraph of this paper one will therefore 
need to eliminate effects due to encounters in the in- 
teraction range which contributes to the TV-invariant 
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c-folding time, but do not, for sufficiently large- 
TV, significantly change particle trajectories on such 
timescales. When this is done, the mathematical 
conditions for large iV-body systems to be described 
by the CBE, for times much longer than the dy- 
namical time, are then also satisfied. An important 
consequence therefore is, if these large-iV systems 
are adequately described by steady state solutions of 
the CBE, the divergence timescale in systems where 
the smoothed out background potential is separa- 
ble will have to increase with N. Moreover, since 
in such a description there are no particle-particle 
correlations and trajectories are independent, diver- 
gence timescales in the full iV-body problem and 
their variation with N should be similar to those of 
trajectories in compatible systems of non-interacting 
background particles. It is the object of this paper 
to examine the validity of such assertions. 

By introducing a short range cutoff, GHH stud- 
ied the behaviour of the divergence timescale due 
to random long (compared to R/y/N) range inter- 
actions. They found it to increase as N 1/3 . Cal- 
culations of the Liapunov exponents from full N- 
body simulations of softened particles over short 
timescales confirmed this for low TV (up to a few 
hundred). One of the goals of the present paper is 
the the investigation of the variation of the exponen- 
tiation times of softened iV-body systems to higher 
N and longer times (Section 0). We confine our- 
selves to spherical systems which, being separable 
in the continuum limit, should have, in that limit, 
steady states completely characterised by regular or- 
bits with no exponential divergence. In addition, in 
an attempt to isolate effects that depend on the full 
self-consistent iV-body interaction from those asso- 
ciated with discreteness noise resulting from rapid 
spatial and time variations of the potential, we also 
examine the divergence of single trajectories moving 
in systems of non-interacting softened particles. We 
do this for three different configurations: a collection 
of non-interacting particles moving in a spherical 
box (Section^), a statistically homogeneous distri- 
bution of fixed particles, and a distribution of par- 
ticles with density decreasing with radius as 1/r 2 
(Section ^J). This latter system is known to contain 
chaotic orbits when non-spherical perturbations are 
introduced, even in the continuum limit. We here ex- 
amine whether this feature leads to different stabil- 
ity properties for its trajectories in comparison with 
the homogeneous case (where, in the aforementioned 
limit, the potential is harmonic and all trajectories 
are regular even in the triaxial case), in the hope of 
exploring the effects of the smoothed out mass dis- 
tribution on the stability of trajectories in discrete 
systems. We summarise and attempt to interpret the 
results in Section [B|. 



2 iV-BODY SYSTEMS 

2.1 Numerical evaluation of Liapunov 
exponents 

Liapunov exponents measure the linear stability 
along the trajectories of dynamical systems. For ex- 
ample, let 



x ; = fi 



(2) 



be the Newtonian equations of motions of an N- 
body gravitational system (with i — 1,N), and 



<5Xi = <5fi 



(3) 



the corresponding variational (i.e., linearised) equa- 
tions. The latter measure the deviation of nearby 
states to the one determined by the first set of equa- 
tions. Next define 



X 



xi, 



, Xn, Xi, 



(4) 



as the 6iV-dimensional "vector field" of the system, 
and 



£ = (<5xi, <5x N , Ski, <5x N ) 



(5) 



as the corresponding vector field of variations. Then 
the above equations can be written as 



X = F 

and 

£ = SF. 



(6) 



(7) 



Now if we choose a particular trajec- 
tory X = X(t,t ,X ) against which we would like 
to measure the deviation, with X being the point 
where we start measuring that deviation (i.e., the 
initial conditions), then (Ym can be rewritten as 



£ = D x F(X(t,t ,X ))£, 



(8) 



where D X F is the Jacobian 6N x 6N matrix 
aFi/Sxj and i,j = 1, 6iV. Now let 



X s =X.(X(i 1 i ,Xo)) 



(9) 



be the fundamental solution of this matrix with the 
initial condition being the identity matrix, the solu- 
tion of (^) is then given by (Wiggins 1991) 

I = x.(t)€o, (io) 

which describes the evolution under the linearised 
dynamics with initial conditions £o in the space of 
linear variations. 

A Liapunov exponent is the infinite limit of 
the "time dependent Liapunov exponent" (Wiggins 
1991) at Xo in the direction £o at time t which is 
given by 



A(e,t) 



€(*) 



■ log 



X.(t)£ 



Coll * °V Ho II 
The Liapunov exponents are then defined as 



(11) 
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ff(£ ,Xo)= lim A(€ ,Xo,t). (12) 

t — >oo 

Numerically of course only A can be calculated. We 
will refer to the inverse of this time dependent Li- 
apunov exponent as the "exponentiation time", the 
"exponential timescale" or the e-folding time. 

For a Hamiltonian system with / degrees of free- 
dom, there are 2/ linearly independent directions in 
phase space for the vector £o to point at, hence there 
are 2/ Liapunov exponents. A positive Liapunov ex- 
ponent indicates unstable behaviour characteristic 
of chaotic motion. Thus, determining the maximal 
exponent is sufficient for detecting the presence of 
such behaviour. The evaluation of the maximal ex- 
ponent is straightforward enough. This is because 
exponential instability, if it is present, will cause al- 
most all initial linear tangent space vectors to realign 
themselves along the subspace of maximal expan- 
sion. A numerical determination of a Liapunov expo- 
nent from almost any initial chosen direction for the 
linear variations will thus tend to give an evaluation 
of the maximal exponent (Wolf et al. 1985) . The only 
complication that arises is that, when the exponen- 
tially increasing solutions of the linearised equations 
become too large, the calculation is slowed down 
(eventually leading to numerical overflow). This is 
easily remedied however by application of the "stan- 
dard algorithm" of Benettin et al. (1976). This algo- 
rithm is based on the local averaging of the deviation 
between neighbouring states, which is done by divid- 
ing the time we run the system into n subintervals. 
An initial linearised deviation £o will therefore be 
transformed into 

-X a £o ? -X s Al s £o ■ ■ ■ ■ ■ -X a . • a J£ g g i£o at times 
t\, t2, t n . The right hand side of (|l l| ) can then 
be rewritten as: 

|| X?...X?XiXX || / || £o ||. (13) 
We now successively define 

& = XUi-i =11 6-i II (14) 
with 

i_! =«*-!/ || 6-1 || (15) 

this means that 

ii x n ...x 2 x X £o n= nr=i II & II ( 16 ) 

and therefore 

t/At 

^standard = hm > . (17) 

n — >og ' ** t 

1 

In practice, this procedure consists of renormalising 
the linearised vector to unity at intervals At, adding 
the logarithm of its norm to the pre-existing sum 
and restarting the integration with this renormalized 
unit vector serving as initial condition for the varia- 
tional (linearised) equations. This avoids numerical 
blowup. 



2.2 Parameters and initial conditions 

In general, the values of the Liapunov exponents will 
depend on the initial conditions. If the phase space 
is "topologically transitive" however — that is if tra- 
jectories can visit any region of the phase space — 
for almost all initial conditions in the connected re- 
gion all the Liapunov exponents will be equal. This 
however is not true of the time dependent Liapunov 
exponents which can be different at any given time. 
Convergence indicates that an invariant phase space 
distribution has been reached. The convergence time 
is therefore a measure of the transport properties of 
a given system's phase space. The time dependent 
Liapunov exponents of systems with more complex 
phase space take longer to converge. Liapunov ex- 
ponents of regular orbits converge to zero at a rate 

lnt 
~ i ' 

Strictly speaking, Liapunov exponents of open 
gravitational systems do not converge at all, because 
these systems are thermodynamically unstable and 
continuously evolve towards more and more con- 
centrated states (e.g., Antonov 1962; Lynden Bell 
& Wood 1968; Padmanabhan 1990). Moreover, the 
evolution rate will depend on N, making the com- 
parison of systems with different iV ambiguous. This 
is not the case however for enclosed A^-body systems 
of single mass particles started from appropriate ini- 
tial conditions and whose energy, mass and radius 
satisfy (El-Zant 1998) 

ER/GM 2 > -0.335. (18) 

We define our system parameters such that 
G — M — R — 1, randomly distributing our parti- 
cles inside R = 1. The quantity on the left hand side 
of condition flis| ) is then solely determined by the 
initial velocities. These are taken to be random in di- 
rection and their absolute values vary with radius as 
v = vo exp(r), where Vo is determined by the initial 
kinetic energy of the system. Since we have already 
determined the spatial distribution, this amounts to 
fixing the virial ratio. Our systems must have virial 
ratios greater than 0.5 if they are to satisfy rela- 
tion ( p_8| ) and be thermodynamically stable. At the 
same time the virial ratio should not be too large, 
so as not to modify the dynamics significantly. We 
choose a virial ratio of 0.69. Systems starting with 
these initial conditions were shown (El-Zant 1998) to 
quickly evolve towards quasisteady isothermal dis- 
tributions, with the total change in the virial ra- 
tio throughout the evolution being of the order of a 
few percent. Systems are enclosed using the method 
also described in El-Zant(1998): particles venturing 
beyond R = 1 are subject to a force of the form 
mK(l - r) 3 , with K = 1300 and m = 1/N is the 
mass of a particle. 

We numerically integrate the Newtonian equa- 
tions of motion along with their variational counter- 
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TIME 

Figure 1. Time variation of inverse of the time dependent Liapunov exponents (calculated from Eq. for different 
particle numbers and random realisations of self gravitating TV-body systems. 



parts starting with random initial conditions for the 
latters (each component of the vectors in the tangent 
space of variations takes a value randomly chosen be- 
tween zero and one and the vector is normalised so 
that its norm is one). In the simulations presented 
here the boundary force is not included in the varia- 
tional equations, however its inclusion was not found 
to significantly affect the results. The integration is 
advanced using a highly accurate variable order vari- 
able stepsize Adams method, as implemented in the 
NAG routine D02CJF with tolerance 10~ 5 . For sys- 
tems consisting of a 100 particles, the dependence 
of the macroscopic evolution on the tolerance was 
checked and found not to vary significantly for toler- 
ances of 10 -3 to 10 -10 . For the tolerance used here, 
the energy was conserved to better than one part 
in 10000 for a hundred units and better than one 
percent for 20000. The energy conservation is much 
better for larger particle numbers. The renormali- 
sation time for the variational equations was taken 
to be one unit. For the above parameters this is of 
the order of a dynamical time. More precisely, if one 
defines the crossing time of a system in virial equi- 
librium by 

T cross = M 5/2 /(2£0 3/2 , (19) 



it is found to be about 4 time units for systems with 
the above parameters. In practice because our sys- 
tems have a virial ratio of 0.69 instead of 0.5, the 
crossing time is 0.67 times shorter. In all simula- 
tions the (Plummer) softening of length is fixed at 
0.1 (> R/y/N for all TV investigated). 



2.3 Results 

For the full TV-body problem, a total of 13 runs were 
investigated. Two different realisations for each TV = 
100, 200, 400, 800 system, three for TV = 1600 and 
four with systems of TV = 3200 particles. The TV = 
100 runs were integrated up to 20000 time units and 
so was one of the TV = 200 runs, the other one was 
stopped at TV = 6600 when it was clear the results 
were converging. The TV = 400 runs were integrated 
up to 7000 time units and the TV = 800 for 5000 and 
10000 units. For TV = 1600, two runs were pursued 
until 4400 units while one was stopped at 1500 units. 
Finally, two of the TV = 3200 runs were integrated 
up to 700 units and the other two for an additional 
400 units, for a total of 1100 time units. 

As can be seen from Fig. [j], for simulations 
that were run up to 4400 time units, the exponen- 
tial divergence timescales show good convergence, 
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400 500 600 700 800 900 1000 1100 

TIME 

Figure 2. Time variation of the inverse of the time dependent Liapunov exponents for two different realisations of 
TV = 3200 self gravitating systems. After the initial evolution, leading to a more or less quasi-steady state, weakly 
damped modes can remain for at least a few hundred dynamical times. 




4 5 6 7 8 



LOG N 

Figure 3. Variation of the exponential timescales with 
particle number for self gravitating TV-body systems. The 
plotted line, which fits well the results for runs with TV = 
100, 200 and TV = 400 has slope of 0.4 



although there seems to be a small systematic de- 
crease in some of the higher-TV runs. This appears 
to be a result of slight increase in central concentra- 
tion accompanying slow evolution towards the in- 
variant equilibrium state (reached much quicker in 
the lower- TV runs). Also, in these runs, it is possible 
to detect small regular oscillations in the exponential 
timescales for early times (smaller than 800 units). 
These are eventually damped out (although varia- 
tions on longer timescales persist for much longer). 
For TV = 3200 however, these oscillations are much 
larger and are not completely damped over the 
timescales investigated (Fig. H). These oscillations 
reflect weakly damped global density oscillations. By 
inspection, at least three timescales can be detected. 
They are of the order of one, ten and few tens of time 
units. 

Already from from Fig. |l| and Fig. ^, it is appar- 
ent that the exponentiation timescale increases with 
TV, but the rate of this increase decreases for larger 
TV. This is evident in Fig. [| where we have plotted 
the final exponentiation times for all the runs as a 
function of TV. The straight line in that figure fits the 
first three values of TV. It has a slope of 0.4, which 
is roughly compatible with a slope of a third, pre- 
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Figure 4. Time variation of the inverse of the time dependent Liapunov exponents (calculated from Eq. for systems 
where one (test) particle interacts with all the other (background) particles, which are non-interacting. The upper panel 
shows results for random initial realisations of systems consisting of 1000 and 2000 particles, while the lower panel shows 
the corresponding results when ./V = 4000, 8000. Invariably, the final exponentiation times are longer for systems with 
larger N. Convergence however is much slower for larger systems (note the different time scales in the upper and lower 
plots). 



dieted by GHH on the basis of their analysis of the 
divergence in a infinite homogeneous medium of soft- 
ened particles, and confirmed by their simulations at 
lower N. One therefore naturally suspects that the 
discrepancy results from self gravitating modes lead- 
ing to fluctuations in the mean field, as evidenced by 
the oscillation of the values of the exponential times 
for higher values of N noted above. These global 
oscillations appear to be quickly washed out by dis- 
creteness noise in the low N (< 800) limit. 

For purposes of comparison, in order to investi- 
gate further the possibility that effects due to self- 
gravity are indeed responsible for the saturation of 
the exponentiation times, we consider in the next 
sections several systems where self gravity is not in- 
cluded. 



3 NON-INTERACTING MOVING 
PARTICLES 

In this approximation one particle interacts with 
all other particles, which do not interact with each 
other (but the motion of which is affected by the 
interaction with that particle). The system is en- 
closed as in the previous section and the total ki- 
netic energy is scaled so as to equal that of the self- 
gravitating systems described there. The interacting 
particle is picked randomly from the spatially homo- 
geneous initial particle distribution, and the full set 
of 12N first order differential equations, describing 
the motion of all N particles and its variations, is in- 
tegrated. The Liapunov exponents are then obtained 
in the same manner as in the previous section (with 
the force law now only containing the interactions of 
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Figure 5. Variation of the exponential timescales with 
particle number for systems similar to those described in 
Fig. ^ (see also text). Least square fits result in power 
laws with index 0.64. 



one particle with all others). All system parameters 
and units (including the softening) are the same as 
before. 

We have run simu- 

lations with TV = 100, 200, 400, 800, 1000, 2000, 4000 
and 8000. For TV = 100, 200 and 400 respectively, 
we integrated 100, 88 and 38 initially statistically 
identical systems for 1100 time units. For TV = 800 
and TV = 1000 we integrated 16 such systems for the 
same time interval. For TV = 2000 there were five 
runs that went for times up to 1100, one up to 6600 
and four up to 25000 time units. Five TV = 4000 runs 
were completed, the first three were integrated for 
8400, 8800 and 10000 time units. These were sup- 
plemented by two 20000 time unit runs. Finally, for 
TV = 8000, there were also five runs, three for 4000 
units and two for 20000. 

Since the field particles are moving, the test par- 
ticle's energy is not conserved. In principle therefore, 
the whole allowed phase space is shared by any two 
initial conditions, provided that the motion is er- 
godic. The Liapunov exponents (and the associated 
exponentiation times) should therefore converge to 
the same asymptotic values for all initial conditions. 
Because the background particles in this system are 
non-interacting, no global modes can be present and 
the system must be in a statistically steady state 
(characterised by a homogeneous density distribu- 
tion). Evolution is then driven by discreteness noise 
alone, the power of which decreases with TV. Diffu- 
sion timescales, therefore, should be longer for larger 
TV, and the convergence timescales should then also 
scale accordingly — since longer time is needed to 



reach an invariant phase space distribution. This 
was indeed found to be the case, as is illustrated in 
Fig. ^ where we show the inverse of the time depen- 
dent Liapunov exponents for pairs of runs with dif- 
ferent particle numbers (1000, 2000, 4000 and 8000). 
Note that both the convergence and exponentiation 
timescales are much shorter for smaller TV. This is 
what is to be expected if the divergence timescale is 
related to "mixing" and the loss of memory of ini- 
tial conditions. Nevertheless, for the runs described 
in this section, convergence was always reached — 
which implies that all trajectories are chaotic and, 
for a given TV, share the same region of phase space. 
As we will see in the next section, this does not ap- 
pear to be always the case with trajectories of non- 
interacting softened systems. 

The variation of the final exponential timescales 
as a function of particle number, for all the runs, is 
shown in Fig. |^. As can be seen, not only is there 
no clear sign of saturation in the increase of the ex- 
ponentiation timescales for the larger TV runs, but 
rate of increase in the exponentiation times is ob- 
viously very different from the self gravitating case. 
The observed TV- variation is very well fit by a power 
law. Best least square fits gave a slope of 0.64 (plus 
or minus 0.008). This is rather different from the 
slope of ~ 1 expected if there was direct correspon- 
dence with the TV variation of the relaxation time 
inferred from two body relaxation theory — even 
though the situation here is clearly similar to the 
idealised systems that are the subject of this theory 
(see, e.g., Chandrasekhar 1942). It is however much 
larger than the slope of the corresponding relation 
for the exponentiation timescales in the full TV-body 
problem — even for small TV. 



4 NON-INTERACTING FIXED 
PARTICLES 

Integrating the full equations of motion, even for en- 
closed systems of non-interacting particles, for long 
intervals, is time consuming. To investigate the be- 
haviour of the Liapunov exponents for larger TV 
we therefore revert to an even simpler situation, 
whereas we integrate individual trajectories moving 
in the potential of systems of fixed particles for short 
times but for a large number of initial conditions. 
For each initial condition, only six first order equa- 
tions of motion along with their associated varia- 
tional equations are integrated. The Liapunov expo- 
nents are then obtained by the same method as in 
the previous sections. 

Two types of distributions are discussed: ho- 
mogeneous and p ~ 1/r 2 . All system parameters 
and units are the same as in the previous sections. 
The velocities are chosen such that the total en- 
ergy is smaller than the energy of the zero veloc- 
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Figure 6. Distribution of the values of the inverse of the time dependent Liapunov exponents after 420 time units 
for 500 trajectories moving in the potential of (constant density) fixed particle systems of N = 100 (upper panel) and 
N = 10000. In the latter case, the bin centered around 250 also contains all trajectories with larger exponentiation 
times. 



ity surface Vzvs for a corresponding smooth distri- 
bution at R = 1 (which, for large N, roughly co- 
incides with the zero velocity surface correspond- 
ing to the actual particle distribution). That is, 
ea ch velocity com ponent i is determined from Vi — 
\J —2(V + Vzvs) x (1 — 2X), where X is a random 
number uniformly distributed between and 1 and 
V is the particle's initial potential energy as deter- 
mined by its initial spatial coordinates. The (Plum- 
mer) softening length is again fixed at 0.1. 

4.1 Homogeneously distributed systems 

Fig. |i] shows histograms displaying the distribution 
of the exponentiation timescales of samples of 500 
trajectories integrated in homogeneous systems of 
100 and 10000 softened particles, with the initial 
conditions described above and for 420 time units. 



As can be seen, the spread in exponentiation times 
increases significantly for larger TV. No obvious cor- 
relation was found between the values of the Lia- 
punov exponents for a given N and the initial con- 
ditions of the trajectories. In fact, starting from the 
same initial conditions and randomly varying the 
positions of the background particles gave similar 
results. This is indicative of a very complex phase 
space structure that may well be worth studying in 
more detail in the future. Perhaps more important, 
it may also indicate that a transition from a highly 
chaotic regime with a single Liapunov exponent to 
a mixed phase space is taking place. 

In the case of systems without discreteness 
noise, that is when the potential does not contain 
rapid spatial variations as is the case here, mixed 
phase spaces have sets of initial conditions from 
which orbits have zero Liapunov exponents (in the 
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Figure 7. Variation of the inverse Liapunov exponents 
of trajectories of homogeneous systems of fixed point par- 
ticles. All trajectories were integrated for 420 time units, 
except in the case N = 100000 where some trajecto- 
ries were also integrated for 4200 units (the exponents 
for these are shown slightly to the right of the 420 unit 
ones). 



infinite time limit). This appears to be the case for 
some of the trajectories of the N = 100000 systems. 
For these, monotonous decrease of the Liapunov ex- 
ponents (roughly with the expected rate of log t/t) 
was detected for times up to 4200 units. Thus it 
appears that some trajectories actually become reg- 
ular for N = 100000, despite the discreteness of the 
system. Nevertheless, most exponents tended to con- 
verge to some non-zero value when the integration 
was continued for longer times. It is possible that 
Arnold diffusion (see, e.g., Lichteberg & Lieberman 
1983), highly ineffective in potentials without high 
frequency variations, is at work here. Since the asso- 
ciated diffusion timescales could be very long, con- 
vergence is expected to be slow. For higher N conver- 
gence is slower — hence the increase in dispersion. 

Fig [j] displays the the changes in the exponen- 
tiation times as a function of TV for subsets of 50 
randomly chosen trajectories of systems with N 
11)0. 500, 1000, 5000, 10000, 100000. As expected, one 
sees an increase of the average value of the exponen- 
tiation times accompanied by an increase in their 
scatter. Even though it appears there is some sat- 
uration in the average values as N increases, this 
may be due to the relatively short integration times 
(which, as is by now clear, can prevent convergence 
of the Liapunov exponents to their final values). For 
example, in the case of N = 100000, we have inte- 
grated 50 randomly chosen initial conditions for ten 
times longer than the other trajectories (integrated 



for the standard time interval of 420 time units). The 
exponentiation timescales for these longer time in- 
tegrations are shown slightly to the right. It is clear 
that they are, on average, larger. Best fits for the 
values of the exponentiation times as a function of 
N are compatible with those found in the previous 
section for the non-interacting moving particles. 

4.2 Systems with density proportional to 

1/r 2 

As was mentioned in the previous subsection, in the 
case of the homogeneous systems studied there, no 
correlation was found between the values of the ex- 
ponentiation times and the initial conditions. As we 
will now see, this is not the case for systems with 
centrally concentrated particle distributions. 

Fig. ^| shows histograms for a sample of 500 
trajectories integrated in particle distributions with 
N — 100 and N = 10000 and average densities vary- 
ing as ~ 1/r 2 . Evidently, the major difference from 
Fig. ^| is that, along with the increase in the av- 
erage value and the scatter of the exponentiation 
times, their distribution becomes bimodal. These 
two groups turn out to be correlated with two dis- 
tinct sets of initial conditions — which in turn can 
be related to stability characteristics of trajectories 
in smooth non-spherical potentials. Recall that in 
such potentials, in the presence of a central mass 
concentration or a central density cusp, elongated 
orbits passing near the centre tend to be chaotic. 
Most box orbits and elongated loops with low initial 
angular momentum (which may also include low an- 
gular momentum orbits in axisymmetric potentials: 
e.g., Caranicolas & Innanen 1991) will fall in this 
category (e.g., Gerhard & Binney 1985; Schwarschild 
1993; Merritt & Fridman 1996; El-Zant & Hassler 
1998). 

Fig. ^ shows the correlation between the ini- 
tial angular momenta of the integrated trajectories 
and the corresponding exponentiation times. As is 
clear, for N = 100, little correlation exists: the dy- 
namics is dominated by the discreteness noise and 
trajectories can diffuse freely from one distribution 
to another. However, for N = 10000 and N = 50000, 
a clear correlation emerges between the two quanti- 
ties, with trajectories starting with higher values of 
the initial angular momentum being more regular. 
Thus, in these cases, the discreteness noise acts as 
a perturbation — having similar effects to those of 
non-spherical perturbations in smooth systems. In 
both situations, trajectories starting at higher an- 
gular momenta are relatively immune to the desta- 
bilising effect of the non-integrable perturbation, as 
compared to the ones starting from low angular mo- 
menta (for an explanation of this phenomenon see 
El-Zant & Hassler 1998). 

In a smooth spherical system, the angular mo- 
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Figure 8. Distribution of the inverse of the Liapunov exponents after 420 time units for 500 trajectories moving in the 
potential of fixed particle systems with p oc 1/r 2 and with N = 100 (upper panel) and N = 10000 



menta of trajectories are exactly conserved. Along 
with energy conservation, these guarantee the inte- 
grability of steady state systems, in turn ensuring 
that all Liapunov exponents tend to zero. The exis- 
tence of non-zero exponents thus implies that angu- 
lar momenta are not conserved. Moreover, if these 
exponents contain meaningful quantitative informa- 
tion about the degree as to which dynamical proper- 
ties differ from the smooth integrable case, then the 
conservation of the angular momenta should corre- 
late with the exponentiation times of our trajecto- 
ries. In order to verify this, the angular momenta 
were sampled at intervals of one time unit, and the 
standard deviation calculated and divided by the 
mean. The results are shown in Fig [lo| where we 
have plotted the exponentiation times as a function 
of the relative dispersion in the total angular mo- 
menta for N = 10000 and N = 50000. In order to 
examine the effect of sampling and to verify conver- 
gence, for the case of N — 10000, we have calculated 



the dispersions over 420 (circles) and 840 (crosses) 
time units (for N = 50000 all results are shown for 
trajectories integrated up to 420 units). The disper- 
sions were found to be nearly identical in both cases, 
but since the Liapunov exponents take longer time 
to converge to smaller values, there is a correspond- 
ing increase in the e-folding times when the integra- 
tion time was doubled. Evidently, for N = 10000 and 
N — 50000, where the motion is not dominated by 
discreteness noise, a clear correlation exists between 
the dispersion and the e-folding times. 

The two sets of trajectories with markedly dif- 
ferent Liapunov numbers turn out to also have dif- 
ferent variation of these numbers as a function of 
N. This is shown in Fig |ll]. The two lines drawn 
are least square fits with slopes 0.28 and 0.72 - 
roughly consistent with those of A^-body systems 
with relatively small N and trajectories in homoge- 
neous distributions of non-interacting softened par- 
ticles respectively. One can also see that the upper 
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LOG OF MOMENTUM DISPERSION 

Figure 10. Relative dispersion in the angular momenta along trajectories of systems of 10000 (upper panel) and 50000 
particles distributed such that p ex. l/r2. In the upper panels the little circles correspond to quantities averaged over 
420 time units while the crosses correspond to ones averaged over 840 time units. 



line fits the data only for N < 10000, for larger N 
saturation appears to set in. Since this could be due 
to the short standard integration time of 420 units, 
we integrated the low Liapunov exponent trajecto- 
ries in systems of N = 10000, 25000, 50000 for up 
to 4200 time units. The results are shown in the 
right hand panel of Fig [ll| Although for some tra- 
jectories the exponentiation time-scales continue to 
rise so as to be consistent with the growth observed 
for smaller values of N, some exponentiation times 
actually decrease — their values becoming consis- 
tent with those of the population with the smallest 
exponentiation times rather than with their origi- 
nal group! This is because, if one integrates long 
enough, trajectories can diffuse through phase space 
between the two distributions — fluctuations can 
eventually lead to changes in angular momenta even 
for the more regular orbits. When this decreases, 
relatively regular trajectories can diffuse into the 



highly chaotic region. On the other hand, it was 
found that for some of the trajectories, the time de- 
pendent Liapunov exponents continue to decrease as 
the integration proceeds (up to 4200 units) these are 
therefore candidates for being regular trajectories. 
Whether longer integrations would confirm this or 
whether they too would eventually escape from the 
regular regions (through Arnold diffusion or some 
other mechanism) is an open question. 

Clearly the phase space structure of these sys- 
tems is highly complex and it may also well be worth 
further investigation. It is clear however that for 
short enough timescales, of the order of hundreds 
of crossing time, a certain convergence towards the 
regular behaviour described by the continuum limit 
is apparent. This can be inferred, for example, from 
Fig. where one can see that for N — 50000 rel- 
atively more trajectories have exponentiation times 
approximating that of the most regular orbits than 
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Figure 11. The variation of the exponentiation times with particle numbers for trajectories of systems of fixed particles 
distributed such as p oc 1/r 2 . Lines with slope 0.28 are fits to the exponentiation times of trajectories trajectories with 
lowest initial angular momenta. Lines with slope 0.72 are best fits to trajectories with highest initial angular momenta. 
For each system, out of the of 500 trajectories integrated, those 50 with the highest Liapunov exponents and another 
fifty with the lowest ones are shown. In the panel on the left all trajectories are integrated for 420 time units. On the 
right, trajectories with the largest exponentiation times are integrated for 4200 in the cases of N = 10000, 25000, 50000. 



for TV = 10000. Moreover, this exponential time is 
constant over a large range of initial momenta and 
is essentially fixed by the integration time — during 
that time, the behaviour of the time dependent Lia- 
punov exponents of these trajectories indicates that 
they have properties of regular trajectories (i.e., the 
exponents monotonously decrease with time). 



5 SUMMARY AND CONCLUDING 
REMARKS 

In this paper the stability of the trajectories of soft- 
ened gravitational systems was examined. The ex- 
ponential instability characterising the divergence of 
nearby trajectories in the case of self-gravitating sys- 
tems of point particles is also present in the corre- 
sponding softened systems. However, whereas in the 
case of unsoftened TV-body systems, the associated 
c-folding times do not depend on the number of par- 
ticles, in the case of softened systems, the timescales 
do vary with TV. For TV up to a few hundred this vari- 
ation is a power law N a with a ~ 0.4, roughly con- 
sistent with a — 1/3 predicted by GHH on the basis 
of analysing random encounters of softened parti- 
cles in an infinite distribution. For larger TV (up to 
a few thousand) however, the exponentiation times 
increase much more slowly. This effect appears to 
be connected to large scale time dependent varia- 
tions in the potential that are swamped at small TV 
by discreteness noise. If this is indeed the case, then 



one might expect the situation to be different if self 
gravity is not taken into account . 

To see if the slow rate of increase in the ex- 
ponential times of TV-body systems did indeed arise 
from collective effects, we have integrated trajecto- 
ries in centrally concentrated and homogeneous dis- 
tributions of non-interacting particles, which were 
either moving independently or were kept fixed. In 
the homogeneous case, not only was there no clear 
sign of saturation in the TV-dependence of the e- 
folding times, but the rate of increase was charac- 
terised by a power law with another, much larger, 
value of the exponent (a « 0.7) for TV up to 100000. 
In addition, the dispersion in the values of the Li- 
apunov exponents integrated from different initial 
conditions increases with TV. This is indicative of 
a complex phase space structure, perhaps a mixed 
phase space where regular and chaotic regions co- 
exist. Indeed, for TV = 50000 and TV = 100000, 
what appeared to be completely regular trajectories 
- characterised by time dependent Liapunov ex- 
ponents monotonously decreasing for thousands of 
dynamical times — were found. This phenomenon 
is, of course, completely absent in the case of point 
particles systems (Valluri & Merritt 1999). 

For TV > 10000, trajectories of systems of fixed 
particles that, instead of being homogeneously dis- 
tributed, had radial density profiles varying as ~ 
1/r 2 could be classified into two groups. High an- 
gular momentum trajectories had, in general, larger 
exponentiation times than low angular momentum 
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ones. This situation is analogous to the case of 
smooth centrally concentrated non-spherical poten- 
tials, where eccentric trajectories can be chaotic. It 
therefore appears that, in the present case, discrete- 
ness noise replaces the effect of global non-spherical 
perturbations as the trigger of chaotic behaviour. 

As apposed to the case of point particle systems, 
the divergence timescales of the trajectories corre- 
late rather well with the conservation of the action 
variables (as evidenced by the conservation of the 
angular momenta) along them. That is, as may be 
expected if the exponents tell us something about 
the phase space diffusion, more chaotic trajectories 
(i.e., those with larger Liapunov exponents) conserve 
their action variables less efficiently than the more 
regular orbits (which should approximate ones in the 
smoothed out potential, where these quantities are 
exactly conserved). 

In centrally concentrated systems, the most 
chaotic trajectories and the least chaotic ones 
also have different iV-variations of the exponential 
timescales, with the former variation being a power 
law with a — 0.28 and the latter with a — 0.72. 
These values correspond to those of JV-body sys- 
tems of a few hundred particles and of trajecto- 
ries of homogeneous non-interacting distributions re- 
spectively. This is result is also analogous to what 
is found in smooth non-spherical: in the continuum 
limit, homogeneous systems have harmonic poten- 
tials and contain only regular orbits. The low Li- 
apunov exponent trajectories in the centrally con- 
centrated systems also correspond to orbits that re- 
main regular when non-spherical perturbations are 
introduced. JV-body systems, on the other hand, are 
inhomogeneous. The largest Liapunov exponent in 
an JV-body system might then be expected to cor- 
respond to that of the most chaotic trajectories in 
a corresponding fixed particle system — provided 
there are no collective effects influencing the results. 
This expected correspondence, and its existence in 
small- TV systems, may be taken as further evidence 
that the saturation in the TV-body case might be due 
to such effects: under the action of collective phe- 
nomena the self consistent field is time varying, and 
therefore the saturation does not contradict the pos- 
sibility that such systems obey the CBE — even in 
spherical but time dependent smooth potentials tra- 
jectories can be chaotic. Weinberg (1998) has shown 
that long lived modes can be triggered by discrete- 
ness noise even for particle numbers far larger than 
those of the JV-body systems described here, and can 
lead to evolutionary effects on a timescale shorter 
than that of the unamplified discreteness noise (see 
also Kandrup & Severne 1986). 

From the above discussion it is clear that the lo- 
cal stability of trajectories of softened gravitational 
systems can tell us much more about their dynami- 
cal properties than in the case of systems consisting 



of point particles. Once the effect of the singular- 
ity in the potential is removed, many more inter- 
esting features become apparent, and obvious cor- 
relations with some important dynamical properties 
are then revealed. It thus appears that the investi- 
gation of the local stability of motion in this case 
should provide clues concerning the mechanisms 
driving chaotic behaviour in gravitational systems 
and their interaction. They include high frequency 
discreteness noise, larger scale collective phenomena 
and the effect of the shape of the global smoothed 
out potential. These are precisely the same phe- 
nomena that drive macroscopic evolution in grav- 
itational systems. Since, in the absence of chaotic 
behaviour, orbital action variables are conserved, 
once "phase mixing" in the corresponding angle vari- 
ables has taken place, no further evolution is possi- 
ble. Therefore, as has been argued in El-Zant(1997), 
the macroscopic evolution of gravitational systems is 
necessarily driven by chaotic motion resulting from 
the aforementioned mechanisms and their interac- 
tions. The study of the local stability of trajectories 
of gravitational systems is thus not just interesting 
on intrinsic merit, but helps isolate and quantify evo- 
lutionary mechanisms. 

Since actual astronomical systems interact 
through the singular Newtonian potential and not a 
softened version of it, one may wonder as to the ap- 
plicability of the results described here to the orbital 
stability of real systems. Nevertheless, most studies 
of large dissipationless stellar systems assume (ei- 
ther explicitly or implicitly) that these systems are 
described by the CBE — an assumption that, as 
pointed out in the introduction, can only, strictly 
speaking, be justified in the case of softened sys- 
tems. The study of softened systems therefore is im- 
portant for testing the validity of this widely used 
assumption. In this paper, as far as I am aware, first 
evidence is given of the convergence towards the reg- 
ular dynamics described by the CBE for steady state 
systems with separable potentials. The fact that this 
convergence is much faster in systems where self 
gravity has been artificially suppressed, and a statis- 
tical steady state forced, may be of crucial physical 
relevance. Nevertheless, if one assumes, as argued 
in the introduction, that the mechanism leading to 
the short TV-invariant e-folding time in point particle 
systems is physically unimportant, then the results 
presented here may be seen as an important step in 
resolving the long standing apparent paradox con- 
cerning the exponential divergence of trajectories of 
gravitational systems. 

In El-Zant (1997) as well as in a couple of other 
papers (El-Zant 1998b; El-Zant & Gurzadyan 1998) 
another, less direct (but perhaps more powerful), ge- 
ometric method chracterising the local divergence of 
geodesies on the (Lagrangian) configuration space 
of dynamical systems, was applied to dynamic and 
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static gravitating ones. There too, attempts to find 
relations between macroscopic evolution and micro- 
scopic instability were made, and correlations be- 
tween some characteristic physical quantities (e.g., 
rotation, clustering and central concentration) were 
also found. In the more abstract geometric setting, 
however, it is more difficult to isolate the physical 
mechanisms at work than when applying the much 
more direct method used in this paper. Moreover, 
in the case of point particles, where in the geomet- 
ric case the divergence is mainly due to the neg- 
ative curvature of the configuration manifold, the 
relation between the methods is far from clear. Not 
only does the negative curvature of the configuration 
space of point particle systems measure time inde- 
pendent deviations between orbits, and not tempo- 
ral trajectories as in the method used here, but also 
the terms giving rise to the exponential instability 
(on a timescale similar to that found using Liapunov 
exponents) involve sums of the first derivatives of 
the potential, and not second derivatives as the case 
when integrating the variational equations. In addi- 
tion, highly technical questions on whether the local 
exponential divergence of geodesies on its own jus- 
tifies the usual global conclusions required to infer 
macroscopic evolution (outside the idealised domain 
encountered in mathematical studies where this was 
shown to be the case: see, e.g., Szczesny & Dobrowol- 
ski 1999) are involved. This is especially relevant 
since point mass systems have singular potentials, 
which makes the metric also singular and the curva- 
ture undefined at some points of the configuration 
manifold, which then becomes "incomplete" (e.g., 
Abraham & Marsden 1978). The global structure of 
geodesies can be affected by these singular "bound- 
aries" . When, by softening the potential, the singu- 
larity is removed, the curvature is no longer nega- 
tive. Chaotic behaviour can then result only from 
fluctuations in the curvature along the motion, the 
amplitude of which will decrease with particle num- 
bers. These fluctuations can also be related to the 
Liapunov exponents (see Casetti, Pettini & Cohen 
2000; Latora, Rapisarda & Ruffo 1999 and the ref- 
erences therein). It is probably then in this regime 
that correspondence between the two methods can 
be most easilly examined. 
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